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This review is a kinetic theory study investigating the effects of inelasticity on the structure 
of the non-equilibrium states, in particular on the behavior of the velocity distribution in the high 
energy tails. Starting point is the nonlinear Boltzmann equation for spatially homogeneous sys- 
tems, which supposedly describes the behavior of the velocity distribution function in dissipative 
systems as long as the system remains in the homogeneous cooling state, i.e. on relatively short 
time scales before the clustering and similar instabilities start to create spatial inhomogeneities. 
This is done for the two most common models for dissipative systems, i.e. inelastic hard spheres 
and inelastic Maxwell particles. There is a strong emphasis on the latter models because that is 
the area where most of the interesting new developments occurred. In systems of Maxwell parti- 
cles the collision frequency is independent of the relative velocity of the colliding particles, and 
in hard sphere systems it is linear. We then demonstrate the existence of scaling solutions for the 
velocity distribution function, F{v,t) ~ VQ{t)~'^ f {{v / VQ{t)) , where vq is ther.m.s. velocity. The 
scaling form /(c) shows overpopulation in the high energy tails. In the case of freely cooling 
systems the tails are of algebraic form, /(c) ~ c~'^~"', where the exponent a may or may not 
depend on the degree of inelasticity, and in the case of forced systems the tails are of stretched 
Gaussian type f{v) ~ exp[— /?(?;/uo)''] with 6 < 2. 

1 Introduction 

The interest in granular fluids [1] and gases has led to a great revival in kinetic theory of dis- 
sipative systems [2, 3, 4, 5], in particular in the non-equilibrium steady states of such systems. 
A granular fluid [6] is a collection of small or large macroscopic particles with short range hard 
core repulsions, which lose energy in inelastic collisions, and the system cools without constant 
energy input. 

As energy is not conserved in inelastic collisions Gibbs' equilibrium statistical mechanics is 
not applicable, and non-equilibrium statistical mechanics and kinetic theory for such systems have 
to be developed to describe and understand the wealth of interesting phenomena discovered in 
such systems. The inelasticity is responsible for a lot of new physics, such as clustering and spatial 
heterogeneities [7], inelastic collapse and the development of singularities within a finite time 
[8,9], spontaneous formations of patterns and phase transitions [10], overpopulated non-Gaussian 
high energy tails in distribution functions [11, 12, 13], break down of molecular chaos [7, 14], 
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single peak initial distributions developing into stable two-peak distributions as the inelasticity 
decreases, at least in one-dimensional systems [15, 16]. These phenomena have been studied 
in laboratory experiments [13], by Molecular Dynamics [7, 14, 15] or Monte Carlo simulations 
[17, 18, 15, 19], and by kinetic theory methods (see recent review [20] and references therein). 

The prototypical model for granular fluids or gases is a system of mono-disperse, smooth 
inelastic hard spheres, which lose a fraction of their relative kinetic energy in every colUsion, 
proportional to the degree of inelasticity (1 — a^), where a with < a < 1 is the coefficient 
of restitution. The model is a well-defined microscopic N particle model, which can be studied 
by molecular dynamics, and by kinetic theory. The single particle distribution function can be 
described by the nonhnear Boltzmann equation for inelastic hard spheres [2, 3]. 

The present article presents a review of kinetic theory studies, dealing with the early stages 
of local relaxation of the velocity distribution t\ and we avoid the long time hydrodynamic 
regime where gradients in density and flow fields are important. So, we resttict our study to 
spatially homogeneous states. Without energy supply these systems are freely cooling [21, 11, 
12]. When energy is supplied to the system a source or forcing term is added to the Boltzmann 
equation [22, 14, 12, 18], and the kinetic equation allows steady state solutions, which depend on 
the mode of energy supply [22, 14, 12, 18]. 

A freely evolving inelastic gas or fluid relaxes within a mean free time to a homogeneous cool- 
ing state, where it can be described by a scaling or similarity solution of the Boltzmann equation, 
t) ~ (l/fo(i))°'/(?^/?^o(*))- Such solutions depend on a single scaling variable c = w/wo(t) 
where v^it) is the r.m.s. velocity or instantaneous width of the disttibution. This early evolution 
is comparable to the rapid decay of the distribution function to a Maxwell-Boltzmann distribution 
in an spatially homogenous elastic system. However, in systems of elastic particles similarity 
solutions of the nonlinear homogeneous Boltzmann equation do not control the long time behav- 
ior of F(f , t) [23]. The earher studies [21, 11, 12] of these problems were mainly focussing on 
inelastic hard spheres, which is the proto-typical model for inelastic gases and fluids, and on the 
extremely simplified inelastic BGK-models with a single relaxation time [24]. 

More recently simplified stochastic models have been introduced [25, 26, 27] to tackle the 
nonlinear Boltzmann equation, while keeping the essential physics of the inelastic collisions. 
Unfortunately the microscopic dynamics of these stochastic models is only defined for velocity 
variables, and the models can not be studied in phase space using the A?^— particle methods of 
statistical mechanics and molecular dynamics. 

Nevertheless, the recent studies of inelastic Maxwell models [28, 29, 30, 31, 32, 33, 34, 35, 
36, 37, 38, 39, 40, 41] have greatly advanced our understanding of kinetic theory for inelastic 
systems, as well as the structure of the resulting velocity distributions and the significance of 
scaling solutions, which are exposing the generic features of relaxation both in homogeneous 
cooling states, as well as in driven systems. 

The physical importance of these scaling solutions has been demonstrated by Baldassarri et 
al. [28, 29, 30] with the help of MC simulations of the nonlinear Boltzmann equation for one- and 
two-dimensional systems of inelastic Maxwell particles. They found that the solution, t), af- 
ter a short transient time, could be collapsed on a scaling form VQ'^{t)f{v/vQ{t)) for large classes 
of initial distributions F(v,0) (e.g. uniform or Gaussian). Moreover, in one dimension they found 
a simple exact scaling solution of the Boltzmann equation, which has a heavily overpopulated al- 
gebraic tail ~ when compared to a Gaussian. In two dimensions they have shown that the 
solutions of the initial value problem for regular initial distributions (say, without tails) also ap- 
proach a scaling form with over-populations in the form of algebraic tails, /(c) ~ c~''~" with an 
exponent a{a) that depends on the degree of inelasticity a. 
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In the same period Krapivsky and Ben-Naim [31, 32, 33], and independently the present au- 
thors [34, 35, 37] developed analytic methods to determine the scaling solutions of the nonUnear 
Boltzmann equation for freely evolving Maxwell models, and in particular to calculate the expo- 
nent a{a) of the high energy tails /(c) ~ l/c''+'^, as a function of the coefficient of restitution 
a. Using the methods of Ref. [12] the present authors have also extended the above results to 
inelastic Maxwell models driven by different modes of energy supply [36]. Here the high energy 
tails turned out to be stretched Gaussians, /(c) ~ exp[— /3c*] with < 6 < 2. Very recently these 
studies were also generalized to inelastic soft spheres [37, 19] both for freely evolving as well as 
for driven systems, where the over-populated high energy tails turn out to be stretched Gaussians 
with < 6 < 2 as well. This class of models covers both inelastic hard spheres and inelastic 
Maxwell models as special cases. We will only briefly touch upon these models in the concluding 
section of this article. 

Subsequently there appeared also rigorous mathematical proofs of the approach to these scal- 
ing forms for freely evolving inelastic Maxwell gases [39], for inelastic hard sphere gases driven 
by white noise [40], and for both types of systems driven by different types of thermostats [41]. 

This discovery stimulated a lot of theoretical and numerical research in solutions of the Boltz- 
mann equation, specially for large times and for large velocities, as universal phenomena manifest 
themselves mostly on such scales. Why were the first results all for Maxwell models? Maxwell 
models [23] derive their importance in kinetic theory from the property that the collision fre- 
quency is independent of the relative impact speed g, whereas the collision frequency in general 
depends on the speed at impact. For instance, for hard spheres the collision frequency is propor- 
tional to g. 

This property of Maxwell models simpUfies the structure of the nonlinear kinetic equation. 
For instance, the equations of motion for the moments (t;") can be solved sequentially as an initial 
value problem; the eigenvalues and eigenfunctions of the linearized collision operator can be 
calculated, and the {2d — 1)— dimensional nonUnear collision integral in the Boltzmann equation 
can be reduced to a {d — 1)— dimensional one by means of Fourier transformation [26]. In 
subsequent sections we will take advantage of these properties to determine similarity solutions 
of the nonlinear Boltzmann equation for the d— dimensional IMM, and the moment equations 
enable us to study the approach of initial value solutions F{v, t) to such similarity solutions. 

The plan of the paper is as follows: In Sections 2 and 3 the nonlinear Boltzmann equation 
for a one-dimensional gas of Maxwell particles is solved without and with energy input to ob- 
tain scaling solutions. The high energy tails are respectively power laws or stretched Gaussian 
tails. Section 4 is an intermezzo about an inelastic BGK-kinetic equation without or with energy 
supply, for which the high energy tails are analyzed. After discussing in Section 5 the basics of 
the nonlinear Boltzmann equation for our two fundamental d— dimensional inelastic models, i.e. 
inelastic Maxwell models (IMM) and inelastic hard spheres (IHS), we repeat the above program 
for cZ— dimensional models in Section 6 for free cooling, and in Section 8 for driven Maxwell 
models. In the intermediate Section 7 we study the approach of the velocity distribution function 
F{v, t) towards the scaUng solution, and Section 9 gives some conclusions and perspectives. 

2 Freely cooling one-dimensional gases 
2.1 Nonlinear Boltzmann equation 

Some of the basic features of inelastic systems can be discussed using the Boltzmann equation 
for a simple one-dimensional model [25] with inelastic interactions, possibly driven by Gaussian 
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white noise. Let us denote the isotropic velocity distribution function at time t hy F{v,t) = 
F{\v\,t). Its time evolution is described by the nonlinear Boltzmann equation, 

^^(^' *) - *) = ^(^1^)' (2.1) 

where the diffusion term, proportional to d'^/dv'^, represents the heating effect of Gaussian white 
noise with strength D. The collision operator, I{v\F), consists of two terms: the loss term, that 
accounts for the so-called direct collisions of a particle with velocity v with any other particle, and 
the gain term that accounts for the so-called restituting collision of two particles with pre-coUision 
velocities v** and w**, resulting in the post-collision velocities {v,w): 



dw 



-F(v** ,t)F(w** ,t) - F(w,t)F(v,t) 

a 



i{v\F) = y, 

= -F{v, t) + ^J duF{u, t)F , (2.2) 

where p = 1 — q = (l-|-a)/2 and a is the coefficient of restitution with < a < 1. All velocity 
integrations in (2.2) extend over the interval (— oo, -|-oo). Because F is normaUzed to unity, the 

loss term is simply equal to —F(v,t). Equation (2.2) is the nonlinear Boltzmann for the Inelastic 
Maxwell Model (IMM) in one dimension, introduced by Ben-Naim and Krapivsky [25]. 

In a direct collision two particles with velocities v and w collide, resulting in post-collision 
velocities, v* and w*, given by: 

V* = v{a) = ^{1 — a)v + ^{1 + a)w or v* = qv + pw 
w* = w{a) = ^{1 + a)v + ^{l — a)w or w*=pv + qw, (2.3) 

where total momentum is conserved. The restituting collisions are events in which two particles 
with pre-coUisional velocities v** and w** collide such that one of the post-coUisional velocities 
is equal to v. The velocities are given by the inverse transformation of (2.3), i.e. v** = v{l/a) 
and w** = w{l/a). 

The velocity distribution function is normalized such that mass and mean energy or granular 
temperature are, 

J dvF{v,t) = 1; = J dvv'^F{v,t) = yvl{t). (2.4) 

For later reference the normalization of the mean kinetic energy is written down for d— dimensional 
systems, where vo{t) is referred as the r.m.s. velocity, and T = as the granular temperature. 
The evolution equation of the energy [v^) is obtained by multiplying (2.1) with / dvv^, which 
yields, 

j^vl = 4D- 2pqvl (2.5) 

It describes the approach of VQ{t) to a non-equilibrium steady state (NESS) with r.m.s. velocity 
vo{oo) = -^20 /pq. Here the heating rate D due to the random forces is balanced by the loss 
rate, oc pqvg, caused by the inelastic collisions. The energy balance equation (2.5) has an stable 
attractive fixed point solution fo(oo) in the sense that any vo{t), with an initial value different 
from vo{oo), approaches this fixed point at an exponential rate. This one-dimensional model 
Boltzmann equation (2.1) will be extended to general dimensionality d in later sections. 
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An inelastic fluid without energy input, in the so-called homogeneous cooling state, will cool 
down due to the coUisional dissipation. In experimental studies of granular fluids energy has 
to be supplied at a constant rate to keep the system in a non-equilibrium steady state, while in 
analytic, numerical and simulation work freely cooling systems can be studied directly. Without 
energy input the velocity distribution F{v, t) will approach a Dirac delta function 5{v) as t ^ oo, 
and all moments approach zero, including the width vo{t). However, an interesting structure is 
revealed when velocities, c = v/vo{t), are measured in units of the instantaneous width vo{t), 
and the long time limit is taken while keeping c constant, the so-called scaling limit. Scaling or 
similarity solutions of the Boltzmann equation are special solutions of a single scaling argument 
c = v/vo{t), where the normalization of mass requires, 

F{v,t) = {vo{t))-''f{v/vo{t)). (2.6) 

Consequently /(c) satisfies the normalizations, 

J dcf{c) = 1; J dccV(c) = id. (2.7) 

To understand the physical processes involved we first discuss in a qualitative way the relevant 
limiting cases. Without the heating term {D = 0), equation (2.1) reduces to the freely cooling 
IMM. If one takes in addition the elastic hmit g — > 0, the collision laws reduce in the one- 
dimensional case to V* = w,w* = v, i.e. an exchange of particle labels, the collision term 
vanishes identically, so that every F{v,t) = F{v) is a solution, and there is no randomization 
or relaxation of the velocity distribution through collisions, and the model becomes trivial at the 
Boltzmann level of description. However, the distribution function in the presence of infinitesimal 
dissipation {a ^ 1) approaches a Maxwellian. If we turn on the noise {D ^ 0) at vanishing 
dissipation {q = 0), the r.m.s. velocity in (2.5), fo(t) = Vq{0) + 4Dt, is increasing linearly with 
time. With stochastic heating and dissipation (even in infinitesimal amounts) the system reaches 
a NESS, through the balance of energy input and dissipation. 



2.2 Scaling solutions 

After this qualitative introduction into the physical processes, we will study scaling solutions (2.6) 
of the nonlinear Boltzmann equation in a freely evolving system without energy input {D = 0). 
We first take the Fourier transform of the Boltzmann equation (2.1). Because of the convolution 
structure of the nonlinear collision term (2.2) it reduces to 

t) = t) + t)^{kp, t). (2.8) 

This equation possesses an interesting scaling or similarity solution of the form, ^{k,t) = 
(j){kvo{t)), which is the equivalent of (2.6) in Fourier representation. Substitution of this ansatz 
into (2.8) yields then, 

-pqk(l)'{k) = -(l){k) + <l){qk)<l){pk), (2.9) 

where the energy balance equation (2.5) has been used to eliminate dvn/dt. By combining solu- 
tions of the form (1 — sk)e'^^ for positive and negative k, one can construct a special solution of 
(2.9) as: (l){k) = (1 -I- exp[— s|/>;|], with an inverse Fourier transform [28] given by 

/(C) = 2/ 212 - (2.10) 
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By choosing s = l/\/2, one obtains the scaling solution that obeys the normalization (c^) = 
^, as imposed by (2.4). The solution above has an algebraic form; it is even in the velocity 
variable, with a power law tail at high energy, and moments of order larger than 2 are 
all divergent. Notice that the solution (2.10) does not depend on the inelasticity, which is an 
exceptional property of the freely evolving inelastic Maxwell model in one dimension, as we 
shall see in later sections. 



2.3 Moment equations and approach to scaling 

An enormous simplification of both elastic and inelastic Maxwell models, as compared to hard 
spheres or other interaction models, is that the infinite set of moment equations can be solved 
sequentially as an initial value problem. This enables us to investigate analytically, at least to some 
extent, how the general solution F{v,t) of the complex nonlinear Boltzmann equation approaches 
the much simpler scaling solution /(c) of a single scaling variable. 

From general considerations we know already that the solution F{v, t) of the inelastic Boltz- 
mann equation approaches a delta function. So, all its moments must vanish in the long time 
Umit, whereas the moments (c") of the scaling form are constant for n = 0, 2, and divergent for 
n > 2 on account of (2.10). 

In our subsequent analysis it is convenient to define the standard moments m„(t) and the 
rescaled moments jU„(t) of the distribution function as, 

runit) = v^finit) = -J dvv^F{v,t). (2.11) 

The evolution equation of the moments can be obtained by multiplying (2.1) by v"' and integrating 
over V, i.e.: 

dtrin 1 f f , , n-r^t ~ 1'^ 



= -Vfln + 



J dvduv''F{u,t)F {^^—^,t^ . (2.12) 



dt pi 

After a simple change of variables, the gain term in this equation transforms into, 



I j dudw{pw + qu)'"'F{u,t)F{w,t) =Y^p''q'^-^mimn-u (2.13) 



n! 



where we have used Newton's binomial formula. We first observe that particle conservation gives 
mo(i) = 1. Combination of (2.12) and (2.13) yields the moment equations, 



dm 

^ + A„m„ = ^ p^q''-^mimn-i, (2.14) 
"-^ 1=2 

where the eigenvalue A„ is given by, 

A„ = l-p"-g" (2.15) 

Next we observe that for an isotropic F{\v\,t) only the even moments are non- vanishing. So, 
(2.14) only involves even values of I and n. This set of equations can be solved recursively for all 
n. For n = 2 the nonlinear term on the right hand side vanishes, and we find, 

m2(0 = m2(0)e-^2t (2.16) 
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with A2 = 2pq. We note that 777-2 (t) = Vq (t) / 4. Similarly we solve the equation for the fourth 
moment, 

7774(0 = [^4(0) + 5mi(0)l exp[-A4i] - ^ml{0) exp[-A2t], (2.17) 

where the equality A4 — 2A2 = —2p^q^ has been used. The coefficients in (2.17) turn out to 
be independent of p and q. Similarly we can show that the dominant long time behavior of the 
higher even moments is given by m,n{t) ~ exp[— A„t], i.e. they decay asymptotically according 
to (2.14) with the nonlinear right hand side set equal to zero. This is the multi-scaling behavior of 
the velocity moments found by Ben-Naim and Krapivsky in [25, 33]. Consequently all moments 
with 77 > vanish d&t ^ 00, consistent with the fact that the distribution F{v, t) S{v) in this 
limit. 

The approach of F{v, t) to a scaling form in d— dimensional inelastic models was first formu- 
lated and conjectured by the present authors based on their analysis of the time evolution of the 
moments (for inelastic hard spheres and inelastic Maxwell models see [35] ; extended to inelastic 
soft spheres [37, 19]), and was subsequently proven in a mathematically rigorous fashion for in- 
elastic Maxwell models and hard spheres [39, 40, 41]. Rather than presenting the mathematical 
proof it is of more interest from the physics point of view to understand how the physically most 
important lower moments of F{v, t) approach their limiting form, and relate these Umits to the 
moments of the scaling solution (2.10). 

We first observe that the assumed large— t behavior (2.6) of F(7;, t) in combination with (2. 1 1) 
implies that /Lt„(t) jin ^or n = 2, 4. . ., where /x„ = (I/77!) / dcc^f{c). So, the moment 
equations for the rescaled moments follow directly from (2.14) and (2.11) to be, 

+ ln^in{t) = J^P'Q'^-'MtW-lit), (2.18) 
"-^ 1=2 

with eigenvalue 

7n = A„-i77A2 (2.19) 

on account of (2.14) and (2.16). For 77 = 2 one obtains that fi2{t) = 1^2 = ^(c^) = i is constant 
for all times, in agreement with the corresponding moment of the scaling form (2.10). Next we 
consider the fourth moment fi4{t), which follows immediately from (2.17) and (2.1 1) and reads, 

fi4{t) = mi{t)/vl{t) = -\^il+ exp(2p2g2i)[^4(0) + (2.20) 

This solution is indeed positive and finite for all finite positive times, and approaches +00 as t 
becomes large. A similar result for the time dependence of /7n(0 is obtained from (2.18) for 
77 = 6, 8, • • •. This behavior is fully consistent with the exact scaling solution (2.10), of which all 
even moments with 77 > 2 are divergent. 

Now, a curious result follows by considering the stationary solutions of (2.18) by setting 
diin{t)/dt = 0. The equations then reduce to a simple recursion relation, 

n-2 

t^n = — J2 P^Q'^'^i^mn-l (2.21) 
1=2 

with initialization, 112 = \- This yields for 77 = 4, 

/^4 = M2 = (2-22) 

74 oZ 
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where we have used the relation 74 = — 2p^gf^, as follows from (2.15). A negative moment of 
a physical (positive) distribution is clearly unphysical. This unphysical behavior continues for 
higher moments, where one finds in a similar manner that the even moments approach a finite 
value, ^2n oc (— )"+^C„, with alternating signs! 

What is happening here? Clearly, the set of solutions {//2n} of the limiting equation (2.21) 
differs from the long time limit /X2„(oo) of the set of solutions H2n{t) of (2.18). The fixed point 
solution {nn, n = 2, 4, . . .} of (2.21) does not represent a stable attracting fixed point for physical 
solutions, but an unstable/repelling one. The physical solutions /in(i) of (2.18) move for large t 
away from the unstable fixed point {i^n} at an exponential rate, given by exp(2p^g^t). 



3 Driven 1-D gases 

In this section we discuss the scaling form of the velocity distribution in the nonequilibrium steady 
state (NESS) for the one-dimensional IMM for different modes of energy supply, and the results 
are compared with those of the proto-typical dissipative fluid of hard spheres in one dimension. In 
general, a NESS can be reached when the energy loss through colUsional dissipation is compen- 
sated by energy supplied externally, as described in (2.5) for the case of external Gaussian white 
noise. This stationary case is described by the Boltzmann equation (2.1) with dtF{v,t) = 0. 
Moreover, at the end of this section all known results for scaling forms in driven inelastic one- 
dimensional systems will be summarized. 

To obtain a description of the NESS, that is independent of the details of the initial state, we 
rescale velocities, c = v/vo{co), in terms of the r.m.s. velocity vo{oo), and express F{v,oo) in 
terms of the scaUng form /(c) introduced in (2.6). The Boltzmann equation (2.1) then takes the 
form, 

-^^§ = -i^'^f('y'^ = '(clf). (3.1) 

where f 0(00) = V^^/PQ- The integral equation for the characteristic function, = / dce~'^^'^ f{c), 
follows by taking the Fourier transform of the above equation, 

(1 + i pqk^)ct^{k) = <P{pk)<P{qk), (3.2) 

where the nonlinear collision term has been obtained as in (2.8). An exact scaling solution of 
this equation in the form of an infinite product has been obtained by Ben-Naim and Krapivsky 
[25, 33], and by Nienhuis and van der Hart [42]. Here we construct the solution following Santos 
and Ernst in [43]. We multiply (3.2) on both sides with 00 (A;) = 1/(1 + ^pqk^), and write the 
equation as an iteration scheme, 

0n+l = Mk)(l>n{pk)(l)n{qk). (3.3) 

The solution can be found iteratively by starting from (j)o{k), and inserting 4>nik) on the right 
hand side to obtain (k) on the left hand side. The result is the infinite product, 



'DO rn 

^w=nn(i+wc.)"^'"', (3.4) 



m=0 e=o 

m 



where = |^ ^ J and kjni = ap ^"^ with a = ^/2/pq. Thus (f){k) has infinitely 
many poles at A; = ±ikme on the imaginary axis with multiplicity (for a / 0). The velocity 
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distribution, 

1 1"^ 

f{c) = — dke^^'^m. (3.5) 
^vr J-oo 

can then be obtained by contour integration. As /(c) is an even function, we only need to evaluate 
the integral in (3.5) for c > 0. The replacement c ^ |c| gives then the result for all c. By closing 
the contour by a semi-circle through the infinite upper half-plane and applying the residue theorem 
we obtain the exact solution in the form of an infinite sum over poles. This representation of /(c) 
is very well suited to discuss the high energy tail. 

The dominant terms for large |c| correspond the poles closest to the real axis, i.e. to the 
smallest values of kmi- In case a 7^ the two smallest values are fcoo = o and k\i = a /p. 
Consequently, the leading and first subleading term are 

/(c) ^ Aoe-"l"l + ylie-("/f)l"l + • • • , (3.6) 



where a = y2/pq and the coefficients are found as 

\f I \ p2(n+l) ^ ^2(n+l) 



^0 = n <5^p 

A, = 




p 



,2(n+l) _ o2(n+l) 



p2(n+l) ^ g2(n+l) 



1 _ p2{n+\) _ g2(n+l) 



(3.7) 



In conclusion, the scaling form /(c) in driven one-dimensional IMM systems shows again an 
overpopulated high energy tail ~ exp(— a|c|), when compared to a Gaussian. However this 
overpopulation is very much smaller than in the freely cooling case, where /(c) ~ 1/c^ for large 
c. 

Figure 1 compares the asymptotic form /(c) k, ^ge""''^' with the function /(c) obtained by 
numerically inverting ip^k) in (3.2) for a = and a = 0.5. Similar results have been obtained by 
Nienhuis and van der Hart [42] and by Antal et al. [44]. We observe that the asymptotic behavior 
is reached for a\c\ <i 4 if a = and for a\c\ <i 8 if a = 0.5. As a = ^/Ijpq this corresponds to 
velocities far above the r.m.s. velocity. 

Interesting limiting behavior is also found in the quasi-elastic limit. This limit is much more 
delicate and requires some care. If we first take a ^ 1 at fixed |c| and next |c| — > 00 (denoted as 
order A in Table I), the high energy tail has a MaxwelUan form. On the other hand, if the Umits 
are taken in the reverse order, i.e. first |c| — > 00 at fixed a < 1 and then a ^ 1 (denoted as 
order B in Table I), the asymptotic high energy tail is exponential. The crossover between both 
limiting behaviors is roughly described by the coupled limit c — > 00 and g — > with the scaling 
variable w = \c\^ = fixed with q = ^(1 — a) <C 1, and occurs at w; 2± 1. If ^u; < 1 the 
distribution function is essentially a MaxwelUan, while the true exponential high energy tail is 
reached ifw^l. 

The results for the scaling form in the quasi-elastic limit not only depend sensitively on the 
order in which both Umits are taken. They also depend strongly on the coIUsional interaction, i.e. 
on the energy dependence of the collision frequency, as well as on the mode of energy supply to 
the system. 

To illustrate this we have collected in Table I what is known in the literature for the different 
inelastic models in one dimension: (i) hard spheres and (ii) Maxwell models, and for different 
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Figure 1: Logarithmic plot of /(c) versus a\c\ for a = and a = 0.5. The dotted lines are the asymptotic 
forms /(c) « Aoe~°'^'^^ at a = and a = 0.5, obtained from (3.7). 

modes of energy supply: (i) no energy input or free cooling, (ii) energy input or driving through 
Gaussian white noise, represented by the forcing term —Dd'^F{v, t)/dv'^ in Boltzmann equation 
(2.1), and (iii) energy input through a negative friction force oc C^^/|^|' acting in the direction 
of the particle's velocity, but independent of its speed. This forcing term, referred to as gravity 
thermostat in [18], can be represented by a forcing term, ({d/dv){v/\v\)F{v, t), in the Boltzmann 
equation. The results for the gravity thermostat corresponding to order A have been obtained by 
the same method as followed in Ref. [15]. 

Inelastic gases in one dimension exhibit the remarkable property that single peaked initial 
distributions develop into double peaked solutions as the inelasticity decreases [15, 16]. It is 
worthwhile to note that in the quasi-elastic Umit a bimodal distribution, ^[6{c + cq) + d{c — cq)] 
with Co = 1/ \/2, is observed in inelastic hard sphere systems both for free cooling and for driving 
through the gravity thermostat. In inelastic Maxwell models however, this bimodal distribution is 
only observed for the gravity thermostat. 

It is also important to note that in the normalization where velocities are measured in units of 
the r.m.s. velocity, the high energy tail in the driven inelastic Maxwell model is only observable 
for very large velocities, as illustrated in Figure 1 for strong (a 0) and intermediate (a = 
i) inelasticity. In the quasi-elastic limit, where (a 1), the tail is even pushed further out 
towards infinity. This also explains how to reconcile the paradoxical results of exponential large-c 
behavior with the very accurate representation of the distribution function in the thermal range 
(c ~ 1), in the form of a Maxwellian, multiplied by a polynomial expansion in Hermite or Sonine 
polynomials with coefficients related to the cumulants (see Refs. [43]). The validity of these 
polynomial expansions, over a large range of inelasticities with (0 < a < 1) has been observed 
before in [12] for inelastic hard spheres and in [27] for inelastic Maxwell models. Derivations of 
the results, collected in Table I can be found in the original Uterature. 
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State System Order A Order B 



Free cooling 


Hard spheres^' ^ 


i [6{c + Co) + 5{c 


- Co)] e-"!'^! 




Maxwell model^' ^ 




c-4 


White noise 


Hard spheres^' ^ 




e-a|c|3/2 




Maxwell model^' ^ 




g-a|c| 


Gravity thermostat 


Hard spheres^' 


i [5(c + Co) + S{c 
\ [5{c + Co) + 5{c 


- Co)] e— ' 




Maxwell model^' ^ 


- Co)] e-"!'^! 



Table 1: Asymptotic behavior of the 1-D scahng form /(c) in the quasi-elastic hmit for orders A and 
B. The first/second footnote in the second column gives the reference where the result for order A/B was 
obtained. References: ^ Ref. [15], ^ Ref. [12], ^ Refs. [28, 25, 34], * Refs. [12, 15], ^ Ref.[43], ^ Refs. 
[15, 36], ^ Ref.[18], ^ Ref. [36J. 



4 Inelastic BGK-Model 



A brief intermezzo on a very simple inelastic Bhatnagar-Gross-Krook (BGK) model, introduced 
by Brey et al. [24], is presented here to show to what extent this model captures the correct physics 
as described by the nonlinear Boltzmann equation for d— dimensional inelastic Maxwell models. 

A crude scenario for relaxation without energy input has to show that the system cools down 
due to inelastic collisions, and that the velocity distribution F {v,t) approaches a Dirac delta 
function 5 (v) ast ^ oo, while the width or r.m.s. velocity vq (t) of this distribution, defined as 
{v^) = ^dvg, is shrinking. Moreover with a constant supply of energy, the system should reach a 
non-equiUbrium steady state (NESS). These features are implemented in a simple BGK-equation, 
i.e. 

dtF {v, t) - DVlF {v, t) = -LOO [F {v, t) - Fq {v, t)] . (4.1) 

Here uoq is the mean collision frequency, which is chosen here as loq = 1 in order to model 
Maxwell models. The terms —u!oF{v, t) and uoFq{v, t) model respectively the loss and (nonlin- 
ear) gain term. This kinetic equation describes the relaxation of F {v,t) towards a MaxwelUan 
with a width proportional to avo (t), defined by 

Fo{v,t) = {^/Travo) '^exp -(v/avo)'^ = (avo)"'' (c/a) , (4.2) 

where c = v/vo. The energy balance equation follows from (4.1) as, 

dv^/dt = 4D-{1- a'^)v^ = AD - 2jv^, (4.3) 

where (1 — a^) measures the inelasticity. In the case of free cooling {D = 0) the solution is 
vo{t) = fo(0) exp[— 7t]. By inserting the scahng ansatz (2.6) into (4.1) we obtain in the case of 
free cooling an ordinary differential equation, that can be solved exactly [24, 37, 19]. Its high 
energy tail is, 

/(c) ~ A/c'*+'^ with a = 1/7 = 2/(1 - a^), (4.4) 

where an explicit expression for the amplitude A has been calculated. 

As we shall see in Section 6, a similar heavily overpopulated tail, / (c) ~ for d > 1, 

has also be found in freely cooling d— dimensional Maxwell models with ujq = 1, where the 
exponent a (a) takes for small inelasticity the form a ~ I/70 = 4o?/ (l — a^). 
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For the BGK-model driven by external white noise we obtain a NESS with a rescaled dis- 
tribution function F (v, oo) = v^'^ (oo) / {v/vq (oo)) with standard width {(?) = ^d, and the 
rescaled equation for /(c) is obtained from (4.1). Its asymptotic solution for c S> a is then, 



/ (c) ~ exp -Pd 



exp 



-2c/\/l-a2' 



(4.5) 



As we shall see in Section 8, a similar exponential high energy tail is found in the white noise 
driven Maxwell models . 



5 d-Dimensional inelastic Maxwell gases 

In this section we consider the spatially homogeneous Boltzmann equation for inelastic Maxwell 
models in d dimensions without energy input. In most models of inelastic particles total momen- 
tum is conserved in binary collisions, and the models qualify as inelastic fluids. The details of the 
collision dynamics are defined in Figure 2. The nonlinear Boltzmann equation for a d-dimensional 
IMM when driven by Gaussian white noise can again be written in the form, 

^^Fiv, t) - DVlF{v, t) = I{v\F), (5.1) 



where the collision term is. 



I{v\F) = J J dw 



-F(v**,t)F(w**,t) - F(v,t)F(w,t) 

a 



(5.2) 



Here /„(•••) = / dn{- ■ ■) is m angular average over a full d-dimensional unit sphere with 

a surface area = /T{d/2). The factor (1/a) in the gain term of (5.2) originates from the 
Jacobian, dv**dw** = {l/a)dvdw. The direct and restituting collisions are given by: 

V* = V — ^{1 + a){g ■ n)n; w* = w + ^{1 + a){g ■ n)n, 

V** =v-^{l + ^){g-n)n; w** = w + ^{1 + ^){g ■ n)n. (5.3) 

In one dimension the dyadic product nn can be replaced by unity, so that the equations above 
reduce to (2.3). 

To point out the differences between the Boltzmann equation for inelastic Maxwell particles 
and the one for the prototypical inelastic hard spheres we also quote the Boltzmann colUsion term 
for inelastic hard spheres, i.e. 

" 1 



I{v\F) = j j dw\g-n\ 



.F{v** , t)F{w** , t) - F{v, t)F{w, t) 



(5.4) 



Here the collision term contains an extra factor \g ■ n\, and the gain term an extra factor 1/a as 
compared to (5.2). The velocity distributions in these inelastic models with or without energy 
input were recently studied by many authors; in particular the inelastic hard sphere gas in [22, 11, 
12, 17, 18, 15, 36] and the inelastic Maxwell models in [25, 26, 27, 28, 31, 34, 35, 42, 44, 43]. 

We return to the IMM and observe that in any inelastic collision an amount of energy ^ (1 — 
o?)(j^ = pqg^^ is lost. Consequently the average kinetic energy or granular temperature Vq keeps 

decreasing at a rate proportional to the inelasticity j{l — a^) = pq. The balance equation can be 
derived in a similar manner as (2.5), and reads for general dimensionality, 

j^vi = AD-{2pq/d)vl (5.5) 
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Figure 2: Geometry of inelastic collisions, where g = v — w is the relative velocity, the unit vector n 
specifies the point of incidence on a unit sphere around the centre of force, 5 = |g x n| = | smd\ is the 
impact parameter, and 6 the angle of incidence. In inelastic collisions the component g^^ = g ■ n = g cos 6 
is reflected, i.e. ^jj = —ag\\, and reduced in size by a factor a = 1 — 2q = 2p — 1, where a is called the 
coefficient of restitution (0 < a < 1), and a = 1 corresponds to the elastic case. If the total energy in a 
collision is = ^(w^+w^), then the energy loss in an inelastic collision is Ai? = — i(l— a^)(7jj = — pg^jj. 

We will discuss freely cooling systems {D = 0) without energy input, as well as driven systems 
(see Section 7), which can reach a non-equilibrium steady state. In the case of free cooUng the so- 
lution of the Boltzmann equation does not reach thermal equilibrium, described by a MaxwelUan, 
but is approaching a Dirac delta function 5^'^'> {v) for large times. However, the arguments given 
in Section 2 suggest again that t) approaches a simple scaling solution of the form /(c), as 
defined in (2.6) after rescaling the velocities as c = v/vo{oo), with normaUzations chosen as in 
(2.7). 

In the case of free cooling the mean square velocity keeps decreasing at an exponential rate, 
^^o(*) = ^^o(0) exp[— pgi/d], but the distribution function rapidly reaches a (time independent) 
scaUng form /(c), which is determined by the nonUnear integral equation, 

as can be derived by substituting (2.6) in (5.1) and rescaling velocities. 

For freely evolving and driven inelastic hard sphere fluids the scaling solutions /(c) have been 
extensively discussed both in the bulk of the thermal distribution, as well as in the high energy 
tails [17, 12, 37, 19]. 
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6 Scaling in d-dimensional free cooling 



6.1 Fourier transform method 

To determine scaling solutions for free cooling {D = 0), we first consider the Fourier transform 
of the distribution function, = (exp[— ifc • v]), which is the characteristic or generating 

function of the velocity moments (f "), and derive its equation of motion. Because F {v,t) is 
isotropic, $(A;, i) is isotropic as well. 

As an auxiliary step we first Fourier transform the gain term in (5.2), i.e. 

Jdvcxp[-ik ■ v] /gain(«l^) = 
J^J dvdw e-Kp[-ik ■ v*]F{v,t)F{tv,t) = /^^ <I>(/cr?+, t)<I>(A;r?_, t). (6.1) 

The transformation needed to obtain the first equality follows by changing the integration vari- 
ables {v,w) — > {v**,w**) and using the relation dvdw = adv**dw**. Then we use (5.3) to 
write the exponent as A; = fc_ • v + fc+ • tu, where 

fc+ = p{k ■ n)n \k^\ = kp\{k ■ n)\ = k'r]j^{n) 

fc_ =k-k+ = k^Jl-z{k- n)2 = fcr/_(n), (6.2) 

with p= l — q= \{l + a) and z = 1 — q^.'in one dimension 77+ (n) = p and rj- (n) = q, and 
can be replaced by unity. The Fourier transform of (5.2) then becomes the Boltzmann equation 
for the characteristic function, 

dt^k,t) = -^{k,t) + ( ^kr]+{n),t)^ikr]_in),t), (6.3) 

Jn 

where <l'(0,t) = 1 because of the normalization of the distribution function, and the collision 
operator is a (d — 1)— dimensional integral. In one-dimension this equation simplifies to (2.8). 

Next we consider the moment equations. Because F{v, t) is isotropic, only its even moments 
are non-vanishing. If we assume that all moments (v") are finite, then the moment expansion of 
the characteristic function takes the form, 

^k,t) = ^ ^—r^iik ■ vD = J2 i-ikTmnit), (6.4) 

n n 

and ^{k, t) is a regular function of k at the origin. In the equation above the prime indicates that 
n is even, and the moments m„ {t) are defined as, 

mn{t) = f3n{v'')/n\, (6.5) 

where Ps = Ini^ ' for real s is an angular integral over n, given by: 

^ j;/'deisiner-Hcoser ^ r(^)r(|) 
j^/^ deism er^^ r(^)r(i)- 

Moreover, the normahzations (2.4) give mo{t) = 1 and m2{t) = \j32{v'^) = \vQ{t). By inserting 
the moment expansion (6.4) in the Fourier transformed Boltzmann equation (6.3), we obtain the 
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equations of motion for the coupled set of moment equations by equating the coefficients of equal 
powers of k. The result is, 

n-2 

thn + Xnrrin = ^h{l,n-l) mi run-i (n > 2) (6.7) 

1=2 

with coefficients, 

\s = l- h{s,0) - hiO,s) = Ql - vU^) - v'-in)], (6.8) 

where all labels {n, s} take even values only. For later use we calculate A2 explicitly with the 
help of (6.2) with the result A2 = 2pq/d. 

To obtain a scaling solution of (6.3) we set ^{k,t) = (j){kvo{t)) where the r.m.s. velocity is 
obtained from (5.5) with D = 0, and reads VQ{t) = vo{0) e-x.p{—pqt/d) = z;o(0) exp(— |iA2). 
This gives the integral equation for the scaling form (p (k), 

-k>^2k^<l> (k) + (/c) = ikv+) <t> (kv-) , (6.9) 
which reduces for d = 1 to, 

-IX^k^cPik) + cj){k) = cl){pk)cl){qk). (6.10) 
Here <p{k) is the generating function for the moments /in of the scaling form /(c), i.e. 

n ' n 

~ 1- |fe2 + A;V4- A;V6 + ---, (6.11) 

where n is even, /xq = 1 and /X2 = ^P2{c^) = 1/4 on account of the normalizations (2.7) and 
P2 = i/das given in (6.6). 

6.2 Small- A; singularity of the characteristic function 

In the previous section we have obtained the Boltzmann equation (6.3) for the characteristic func- 
tion, the moment equations (6.7), and the integral equation (6.9) for the scaling form (j)(k). These 
equations provide the starting point for explaining the MD-simulations of Baldassarri et al. for 
the freely cooling IMM in two dimensions, i.e. data collapse after a short transient time on a 
scaling form /(c) with a power law tail. In the present section we derive the solution /(c) with a 
power law tail, and in Section 7 we will study the approach to this scaling form. 

The strategy to determine analytically a possible solution with a power law tail is by assuming 
that such solutions exist, then inserting the ansatz /(c) ~ l/c'^+'^ into the scaling equation (6.9), 
and determining the exponent a such that the ansatz is indeed an asymptotic solution. We proceed 
as follows. Suppose that /(c) ~ 1/c""'"'^, then the moments of the scaling form /(c) are 
convergent for n < a and are divergent for n > a. As we are interested in physical solutions 
which can be normalized, and have a finite energy, a possible value of the power law exponent 
must obey a > 2. 

The characteristic function is in fact a very suitable tool for investigating this problem. Sup- 
pose the moment /x„ with n> a diverges, then the n-th derivative of the corresponding generating 
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function also diverges at k = 0, i.e. (p {k) has a singularity at k = 0. Then a simple rescaling 
argument of the inverse Fourier transform shows that 4> (k) has a dominant small-fc singularity of 
the form (j) (k) ~ \k\"', where a / even. On the other hand, when all moments are finite, the 
characteristic function (p{k) is regular at the origin, i.e. can be expanded in powers of k"^. 

We first illustrate our analysis for the one-dimensional case. As the requirement of finite total 
energy imposes the lower bound a > 2 on the exponent, we make the ansatz, consistent with 
(6.1 1), that the dominant small- A; singularity has the form, 

cl){k) = l-lk^ + A\k\'' , (6.12) 

insert this in (6.10), and equate the coefficients of equal powers of k. This yields, 

^aX2 = Xa or apq = 1 — — g". (6.13) 

The equation has two roots, a = 2, 3, of which a = 3 is the one larger than 2. Here A is left 
undetermined. Consequently the one-dimensional scaUng solution has a power law tail, /(c) ~ 
in agreement with the exact solution (2.10). 
For general dimension we proceed in the same way as in the one-dimensional case. We insert 
the ansatz (6.12) into (6.9), and equate the coefficients of equal powers of k. This yields for 
the coefficient of k^ the identity 2pq/d = A2, and for the coefficient of k°- the transcendental 
equation, 

laX2 = JJl-v%-V-]=K- (6.14) 

The equation above obviously has the solution a = 2. We are however interested in the solution 
with a > 2. In the elastic limit (a 1) the solution is simple. There q ^ and a diverges. The 
contributions of 77^ on the right hand side vanish because rf± < 1, and the exponent has the form, 

a^^ = j^. (6.15) 
pq 1 — 

This result has qualitatively the same shape as the numerical solution of (6. 14), shown in Figure 
3 for d = 2. Moreover the simplified BGK-model of Section 4 predicts the same qualitative 
behavior for the exponent of the power law tail. For general values of a one needs to evaluate the 
integrals h{a, 0) and h{0, a), defined in (6.8). 

Here we only quote the results, and refer for technical details to the Uterature [35, 31], i.e. 

h (a, 0) = p'^Pa, h (0, a) = 2F1 (-§, i; |; z) , (6.16) 

where Pa is given in (6.6), and 2F1 is the hypergeometric function with z = I — q'^. One can con- 
veniently use an integral representation of 2F1 to solve this transcendental equation numerically. 
We illustrate the solution method of (6.14) with the graphical construction in Figure 4, where we 
look for intersections of the line y = |sA2 = 70s with the curve y = Xs for different values of a. 

The relevant properties of are: (i) lims_>o A^ = — 1 whereas Aq = because of particle 
conservation; (ii) A,, is a concave function, monotonically increasing with s, and all eigen- 
values for positive integers n are positive (see Figure 4). As can be seen from the graphical 
construction, the transcendental equation (6.14) has two solutions, the trivial one (so = 2) and 
the solution si = a with a > 2. The numerical solutions si(a) for d = 2 are shown in Figure 3 
as a function of a, and the a-dependence of the root a{a) can be understood from the graphical 
construction. In the elastic limit as a j 1 the eigenvalue A2(a) because of energy con- 
servation. In that limit the transcendental equation (6. 14) no longer has a solution with a > 2, 
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Figure 3: Exponent a{a) (solid line) as obtained by numerical solution of (6.14). The open circles repre- 
sent the exponents, measured from the velocity distribution functions obtained from MC simulations, and 
shown in Fig. 5. Figure from Ref.[19]. 

and a(a) oo according to (6.15), as it should be. This is consistent with a Maxwellian tail 
distribution in the elastic case. Krapivsky and Ben-Nairn [31] have also solved the transcendental 
equation asymptotically for large d (d S> 1), which gives quaUtatively the same results as those 
shown in Figure 3 for two dimensions. 

These results establish the existence of scaling solutions / (c) ~ 1/ d^^"- with algebraic tails, 
where the exponent a is the solution of the transcendental equation (6.14). Using a somewhat 
different analysis Krapivsky and Ben-Naim [31] obtained the same results for the algebraic tails 
in freely cooling Maxwell models. 

The previous results have been confirmed in a quantitative manner by means of MC simula- 
tions in [30, 19] for different values of a. The algebraic tails of /(c) are shown in Fig. 5, and 
the exponents a, measured from the MC data in Figure 5 are plotted in Fig. 3. Both graphs show 
excellent agreement with the analytic results, derived here. 

7 Moment equations and approach to scaling 
7.1 JMoments of velocity distributions 

In this section we study the effects on the moments of power law tails in the scaling form, and 
we analyze in what sense the even moments mn{t) = (3n{v"')/n\ at large times are related to the 
moments /U„ = f3n{c^)/n\ of the scaUng form /(c) ~ l/c"+'', which are divergent for n > a 
and xtTanin finite for n < a. 

First consider the moment equations (6.7) where mo (i) = landm2(i) = 777,2(0) exp[—A2i] = 
\vQ{t). Similarly one shows [35, 31] that mn{t) ~ exp[— A„t] for large t. Consequently all mo- 
ments with n > vanish ast —>■ 00, consistent with the fact F{v, t) —>■ S^'^^ {v) in this limit. 
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Figure 4: Graphical solution of (6.14) for different values of the paramenter a. The eigenvalue A., is 
a concave function of s, plotted for different values of the restitution coefficient q for the 2-D inelastic 
Maxwell model. The line y = 570 — is plotted for a = 0.6, 0.8 and a = l(top to bottom). The 

intersections with Xg determine the points sg (filled circles) and si (open circles). Here si = a determines 
the exponent of the power law tail. For the elastic case (a = 1, g = 0) there is only one intersection point. 

The rescaled moments = mn{t)/vQ{t) show a more interesting behavior. We analyze 

how they approach their hmiting form //„(oo). Inserting this definition into the moment equations 
(6.7) and using vo(t) = 'Uo(O) exp[— iA2t] we find for the rescaled even moments with n > 0, 

n-2 

An(*) +7n/^n(i) = ^ h{l,n - l)lJ.i{t)Hn-l{t) 

1=2 

In = K- ^nX2. (7.1) 

The infinite set of moment equations (7.1) for can be solved sequentially for all n as 

an initial value problem. To explain what is happening, it is instructive to consider again the 
graphical solution of the equation, 7^ = — ^3X2 = for different values of the inelasticity 
q or a, as illustrated by the intersections {sq, Si} of the curve y = Xg and the line y = ^3X2, 
where sq and si are denoted respectively by filled (•) and open circles (o). These circles divide 
the spectrum into a stable branch (sq < s < si) and two unstable branches (s < sq) and (s > si). 
The moments //^(t) with s = n > a are on an unstable branch (75 < 0) and will grow for large 
t at an exponential rate, /x„(t) 2± /x„(0) exp[|7„|f], as can be shown by complete induction from 
(7.1) starting at n = [a] + 1. They remain positive and finite for finite time t, but approach +00 
as t 00, in agreement with the predictions of the self consistent method of Section 6. The 
moments /x„ with n = 2, 4, • • • , [a] with n on the stable branch are globally stable and approach 
for f — 00 the Umiting value /tt„(oo) = which are the finite positive moments of the scaUng 
form (7.2), plotted in Fig. 6. In summary, ^ 00 if n > a, and fin{t) approach /x„(oo) = fin 
for n < a, in agreement with the predictions of the power law tails in Section 6.2. 

The behavior of the moments described above is considered as a weak form of convergence or 
approach of the distribution function F{v, t) to the scaling form /(c) for t ^ 00. The physically 
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Figure 5: Simulation of the velocity distribution function by the MC method, showing power law tails. 
Figure from Ref. [19]. 

most relevant distribution functions are those with regular initial conditions, i.e. all moments 
m„(0) = <(0)//n(0) < oo. 

7.2 Moments of scaling forms 

Next we consider the moments jLt„ generated by the scaling form in (6.12), corresponding to 
/(c) ~ l/c"-'^'^, where a > 2 and not equal to an even integer. This implies that the n— th order 
derivatives of ^(A;) at = 0, and equivalently all moments are finite if ra < [a] < a, and all 
those with n > a are divergent. Here [a] is the largest integer less than a where [a] may be an 
even or odd integer. Hence, the small-k behavior of ^(A;) can be represented as , 

[a] 

where the prime on the summation sign indicates that n takes only even values. The remainder 
is of order o(|/c|") as A; ^ 0. In this finite sum we only know the exponent a and the moments 
/xo = 1 and 112 = 1/4. Now we calculate the unknown finite moments of the scaUng form, 
with 2 < n < [a]. This is done by insert ing (7.2) into the kinetic equation (6.9), yielding the 
recursion relation, 

n-2 

IJ-n = (l/7n) XI ^(^' ^ ~ l)l^ll^n-l (7.3) 
1=2 

with initialization ii-2 = 1/4, where / and n are even. The solutions /in for n = 4, 6, 8 are shown 
in Fig. 6 as a function of the coefficient of restitution a. The physical branches of these functions 
are the ones that start positive at a = 1. Furthermore we observe that the root s = a of the 
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Figure 6: Moments ne and /xg as a function of the coefficient of restitution a. Starting at a = 1 
the moment /i„ increases monotonically as a decreases following the physical branch (thick line) untill 
a reaches a zero of 7s, where /z„ diverges towards +00. The recursion relation (7.3) has a second set of 
solutions {/Un} that become negative for small a, indicating unphysical solutions. 

transcendental equation (6.14), 7^ = As — ^3X2 = 0, indicates that 7^ changes sign at si = a (see 
open circles in Fig.4). This change of sign, where the branch becomes unphysical, can according 
to Section 6 be interpreted as all moments with n > si = a {n on unstable branch) becoming 
divergent. That this interpretation is the correct one, has already been demonstrated in the Section 
7.1, where it is shown that forn > a the reduced moments jU„(f) — >^ 00 as t — > 00. 

The recursion relation (7.3) for the moments fin in the one-dimensional case is again a bit 
pathological in the sense that the stable branch (sq = 2 < s < si = 3) contains only one single 
integer label, i.e. s = 3. So only //q = 1 and //2 = 1/4 are finite, and all other moments are 
infinite, in agreement with the exact solution of Baldassarri et al. 

The recursion relation (7.3) has a second set of solutions {fin}, simply defined by iterating 
the recursion relation for n arbitrary large. This set contains negative moments fin [26]. The 
argument is simple. Consider fin in (7.3) with n = [a] + 1. Then the pre-factor l/7[a]_|_i on the 
right hand side of this equation is negative because the label [a] + 1 > a is on the unstable branch 
of the eigenvalue spectrum in Fig.4, while all other factors are positive. This implies that the 
corresponding scaling form /(c) has negative parts, and is therefore physically not acceptable. We 
also note that the moments fin of the physical and the unphysical scaling solution (j){k) coincide 
as long as both are finite and positive in the a— interval that includes a = 1. These unphysical 
solutions are also shown in Fig.6. 
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8 Driving and non-equilibrium steady states 



8.1 Energy balance 

The present section is devoted to the study of systems of inelastic Maxwell particles with energy 
input. Here the system may or may not be able to reach a non-equilibrium steady state (NESS). To 
reach a spatially homogeneous steady state, energy has to be supplied homogeneously in space. 
This may be done by applying an external stochastic force to the particles in the system, or by 
connecting the system to a thermostat, modelled by a frictional force with negative friction. Com- 
plex fluids (e.g. granular) subject to such forces can be described by the microscopic equations 
of motion for the particles, = Vi, and = + (i = 1, 2, • • •), where and are the pos- 
sible friction and random forces respectively. If needed one may also include in a conservative 
(velocity independent) force. 

Regarding the Negative Friction (NF) thermostats, the most important and most common 
choice [18] is a friction, linear in the velocity, +jv, the so-called iso-kinetic or Gaussian ther- 
mostat [45, 46]. A second example of a negative friction force is a = C,v, which is acting in 
the direction of the particle's velocity, but independent of its speed. Furthermore, the external 
stochastic force, ^j, is taken to be Gaussian white noise with zero mean, and variance, 

= 2D5ij5^05{t - t'), (8.1) 

where a, (5 denote Cartesian components, and D is the noise strength. The Boltzmann equation, 
describing a spatially homogeneous system driven in this manner, takes the form, 

dtF{v) + (V^ • a - DVl)F{v) = I{v\F), (8.2) 

where !F = • a — DV^ represents the driving term. 

Next we consider the balance equation for the granular temperature in driven cases, where 
the external input of energy counterbalances the coUisional cooling, and may lead to a NESS. We 
proceed in the same manner as for the free case, and apply (/ dvv^) to the Boltzmann equation 
in (8.2) with the result, 

d{v'^)/dt = J dvv^I{v\F) + 2{va) + 2dD. (8.3) 

The second and third term are obtained from the driving term in (8.2) by performing partial 
integrations. The most common way of driving dissipative fluids in theoretical and MD studies 
[12, 18, 31, 27, 47, 48] is by Gaussian white noise (WN) {a = 0; D 0). We include in our 
studies also the two types of NF-thermostat (a / 0; D = 0), discussed above. The resulting 
energy balance equation is, 

f - >o (WN) 

2CM _ My^ (const NF) (8.4) 
{l-^)vo (iso-kin) 

Here we have used the relation, {\v\) = fo(|c|), where the last average (|c|) is a moment of /(c). 
One sees that the coUisional loss is counterbalanced by the heat, generated by randomly kicking 
the particles or by the negative friction of the Gaussian thermostat. 



dvp 
dt 
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The stationary solutions of the first two equations are stable attractive fixed points, which are 
approached at an exponential rate, 



y/2dD/pq (WN) 
".M = ^ m (const NF) 

The case of driving by an iso-kinetic thermostat (with linear negative friction) is a marginal case, 
discussed in [19], because stationarity is only reached when the friction constant has exactly the 
value 7 = pq/d, in which case any initial value vo{0) is stationary. So, for an inelastic Maxwell 
gas, driven by an iso-kinetic thermostat, there does not exist a stationary state, because the general 
solution of the rate equation for arbitrary value of 7 is, 

vo{t) = vo{0) exp[(7 - pq/d)t\ (iso-kin), (8.6) 

which may be increasing or decreasing as t ^ 00, depending on the inequalities, 7 < 70 or 
7 > 70. 

8.2 Scaling equation 

The equations (8.2) and (8.5) show that the NESS solution F{v, 00) depends strongly on the mode 
of energy supply. To exhibit possible universal features of the solution we measure velocities 
in their typical magnitude, i.e. the r.m.s. velocity wo(oo), just like in thermal equilibrium, and 
introduce the rescaled distribution, 

F{v, 00) = -^f (-^) . (8.7) 
ug(oo) \vo[oo)J 

The integral equation for the scaling form /(c) follows in this case by inserting (8.7) in (8.2), and 
setting dtF = 0, i.e. 

^ :V.V = ~?^ylf (WN) 



I{c\f) = { ' l'^ " (8.8) 

-^"^c-m = ^^^Vc-{cf) (const NF) 

In the second equality on both lines fo(oo) has been eliminated with the help of (8.5). 

Next we consider the special case of a system driven by an iso-kinetic thermostat (a = 
^v;D = 0), where F{vJ.) does not approach a NESS, but rapidly reaches a scaling state, de- 
scribed by (2.6) and having the time dependent r.m.s. VQ{t) in (8.6). In that case the terms dtF 
and TF in Eq. (8.2) produce respectively the terms (—7 + pq/d) Vc • (c/) and 7Vc • (c/). 
So the terms containing the friction constant 7 cancel. The scaling equation for the iso-kinetic 
thermostat then becomes, 

I{c\f) = ^Ve • (c/) (iso-kin). (8.9) 

This scaling equation is identical to the one derived in (5.6) for free cooling, and no trace of the 
friction constant 7 of the iso-kinetic thermostat remains. For the case of inelastic hard spheres the 
equivalence of the integral equations for the scahng form in both cases has been observed before 
by Montanero and Santos [18]. However the big difference between inelastic Maxwell particles 
and inelastic hard spheres is that the latter system has an energy balance equation with a stable 
attracting fixed point solution, but a NESS does not exist in the former case. 
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8.3 High energy tails 

The scaling equations for inelastic Maxwell particles [36] in the previous subsection cannot be 
solved exactly. However its asymptotic solution for c » 1 can be determined by the same proce- 
dure as used for inelastic hard spheres [12]. To do so we neglect the gain term in I{c\f) in (8.8). 
Then /(c|/) is replaced by IiQ^^{c\f) ~ —f{c), and the asymptotic solutions of (8.8) are found 
in the form of stretched Gaussians, /(c) ~ exp[— /3c^] with positive h and [3. One then verifies a 
posteriori that the stretched Gaussians are indeed consistent solutions of (8.8) and (8.9) by sub- 
stitution them back into I{c\f), and showing [36] that the loss term is asymptotically dominant 
over the gain term as long as exponent b and coefficient /? are positive. 

After these preparations we insert the stretched exponential form into the Boltzmann equation 

(8.8) , and match the leading exponents on both sides of the equation, as well as the coefficients in 
the exponents of these terms. This gives the following results for the asymptotic high energy tail, 
/(c) ~ exp[— /3c^], in d— dimensional inelastic Maxwell models, 

h = l /? = vlf (WN) 

6 = 1 /3=^ (const NF) (8.10) 
6 = inconsistent (iso-kin). 

We conclude this subsection about driven inelastic Maxwell models with some comments: 
Why is the result, 6 = 0, inconsistent for the iso-kinetic thermostat in this model? Taking the 
limit 6 ^ in exp[— /3c^] suggests that /(c) has indeed a power law tail. As we have seen in 

(8.9) , the scaling equation for the IMM driven by this thermostat is equivalent with the scaling 
equation for free cooling, and we know from the analysis in Section 6.2 that /(c) ~ l/c"+^ has 
indeed a power law tail. However the exponent a that would have been obtained from (8.9) with / 
replaced by Ij^gg ~ —/(c) does not yield a consistent solution. In the limit 6 ^ the gain term 
^gain ^'^ longer be neglected with respect to the loss term. That this is indeed the case can 
be seen from (6.14), where the terms 77° and rfi originate from the gain term. These terms give 
substantial contributions to the value of a, and are even dominant for small values of a. We also 
note that in the case of driving by an iso-kinetic thermostat — which turns out to be equivalent to 
the free cooling IMM system — ^the driven system does not reach a nonequilibrium steady state, 
but keeps either heating up or cooling down, depending on the thermostat strength 7. 

White noise driving and positive h lead for all d > 1 to consistent asymptotic solutions of 
the scaling equations with overpopulated high energy tails of simple exponential type, /(c) ~ 
exp[— ci/M/m] (Refs.[36, 38]), in complete agreement with the asymptotic result (3.6), and in 
qualitative agreement with the corresponding result (4.5) for BGK-models. Here all moments 
/ dec"- f{c) < 00. This would not be the case for power law tails. The one-dimensional version 
of this problem has been extensively studied by Nienhuis and van der Hart [42], and by Antal et al. 
[44] using MC simulations of the Boltzmann equation. MC simulations of the two-dimensional 
version of this problem have been carried out in Ref. [19]. 

When the IMM system is driven by a constant A^F-thermostat the scaling function also shows 
an exponential tail, /(c) ~ ex.p[—2d{c)c/pq], with a coefficient that depends on the first moment 
(c) of the complete scaling /(c). In Refs. [12, 37, 19] methods have been developed to calculate 
these moments perturbatively. 

For comparison we also quote the high energy tail for d— dimensional inelastic hard spheres, 
which are also of the form of stretched Gaussians, and we quote the results for the exponents b 
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and the coefficients /3 for the different modes of energy supply, i.e. 



^' = 3/2 /3 = ^/if (WN) 

6 = 2 /3 = ^-vfc ('^^"^^NF) (8.11) 
6=1 ^ (iso-kin) 

where /3i is given in (6.6), (c) = / dccf{c) and K3 = /^^ / dcdit|(c — u) • n\^f{c)f{u) [37]. 
For a close to unity the replacement of /(c) by the Gaussian vr"'^/^ exp[— c^] gives a good ap- 
proximation, yielding (c) ~ l/[v^/3i] and K3 = y/^/ir. The results for the WN- and iso-kinetic 
thermostat have been first derived in [12], and confirmed by MC simulations in [18]. The the- 
oretical result for the constant NF thermostat was first derived in [18], but its consistency was 
questioned. For that reason we have verified the consistency of the result (8.11) a posteriori, and 
we confirm that it is indeed an asymptotic solution of (8.8) with the full Boltzmann collision 
operator, as long as a < 1. 



9 Conclusions 

In the present paper we have reviewed the new developments on anomalous velocity distributions 
in gases of inelastic Maxwell models (IMM), and compared the results with those for the proto- 
typical granular model, the inelastic hard spheres (IHS). The velocity distributions, obtained in 
this article, are scahng or similarity solutions (2.6) of the nonlinear Boltzmann equation. 

The dominant common feature in all these results is that the inelasticity of the collisions 
creates over-populations of high energy tails of the velocity distribution F{v, t), when compared 
to the Gaussian Maxwell-Boltzmann distribution F{v,oo) in thermal equilibrium of systems of 
elastic particles. At finite inelasticity (a < 1) the overpopulations in the high energy tails are 
power laws, /(c) ~ l/c""*"^ (free cooling in IMM), or stretched Gaussians, exp[— /3c^] with 

< 6 < 2 (free cooling IHS and driven IMM and IHS). 

An intermezzo, presented in Section 4, shows that the results obtained for inelastic Maxwell 
models are rather robust. We consider an extremely simplified inelastic BGK-model, and show 
that the resulting over-populations of high energy tails, both with and without energy input, are 
qualitatively the same as for the nonlinear Boltzmann equation of inelastic Maxwell models in 
d— dimensions. 

Returning to the main menu of IMM's, we note that the degree of overpopulation is decreasing 
{b t 2) with the increasing efficiency to randomize the velocities of the highly energetic particles, 
either by collisions or by the external forcing terms. As the IHS's have a larger collision frequency, 
oc g, than the IMM's with a constant collision frequency, the IHS have smaller over-populations 
than the IMM's. Because external white noise apphed to freely coohng inelastic gases adds an 
extra mechanism for randomization, the tails in the driven cases show lower over-populations than 
in the freely cooling case. 

A intriguing question about over-populated tails: "Are power laws or stretched exponentials 
the generic form of over-populated tails in inelastic models", is difficult to answer because we 
have only information from two different interaction models. In two very recent articles [37, 
19] new classes of inelastic models have been introduced, that correspond to soft spheres with 
repulsive interactions (1/r^). These soft spheres have a collision frequency oc g'^ with u = 

1 — 2(d — l)/s. The limit s — > 00, or | 1, corresponds to sttong IHS interactions at high energy. 
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In the limit | 0, analogously to s | 2(d — 1), the interactions decrease. For > the IMM- 
interactions {u = 0) are the weakest of all. The results for these inelastic soft spheres models 
[37, 19], corresponding to (8.10) and (8.11), are again stretched Gaussians /(c) ~ exp[— /Jc**] 
with 

6= ^(i^ + 2) (WN; iy>0) 

(9.1) 

b = u (iso-kin; u > 0) 

The limit 6 = | for the iso-kinetic thermostat, which also corresponds to free cooling (see 
Section 8), is consistent with a power law tail. This analysis shows that the generic type of over- 
population is a stretched Gaussian. A freely cooUng IMM with /(c) ~ 1/c""'"'^ is an isolated 
borderline case, that is most likely not the best model to describe the short range, hard core 
impulsive interactions of granular particles. 

At the end of Section 3 we have seen in Table I the phenomenon of compressed Gaussians 
with b > 2, corresponding to under-population of high energy tails. The solutions with two delta 
peaks are extreme cases of compressed Gaussians. More extensive discussions of compressed 
Gaussians and of peak-splitting in one-dimensional inelastic models can be found in [15, 43, 19]. 
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